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Abstract 

We study vortex-lattice phases for a Bose gas trapped in a rotating optical-lattice near the 
superfluid-Mott-insulator transition. We find a series of abrupt structural phase transitions where 
vortices are pinned with their cores only on plaquettes or only on sites. We discuss connections 
between these vortex structures and the Hofstadter-butterfly spectrum of free particles on a rotating 
lattice. 
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I. INTRODUCTION 



Two of the most exciting directions in cold-atom research involve studying lattice sys- 
tems and rotating systems l|. By increasing the importance of interactions compared to 
kinetic energy, lattices allow one to study strongly correlated phenomena such as the boson 
superffuid-Mott-insulator transition 2|. These lattice systems are ideal for studying model 
many-body systems and protocols for quantum information processing [3|. Rotating gases 
lead to interesting vortex physics [4, 0, |y, Q], and the promise of exotic states such as those 
which give rise to analogs of fractional quantum-Hall effects 

aa 

. Here we study the inter- 
play of lattice physics and rotation physics by calculating the vortex-lattice structures near 
a Mott transition. 

In the absence of an optical lattice a rotating Bose-Einstein Condensate (BEC) develops a 
triangular lattice of singly quantized vortices 5|, |6) . This triangular configuration minimizes 
the logarithmic vortex-vortex interaction. However, as seen in recent experiments far from 
the Mott regime 10], a sufficiently deep optical- lattice potential will pin these vortices at 



the maxima of that potential 
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In this paper we show that qualitatively different behavior can be seen in the superfluid 
state near the Mott-insulator phase. We find that due to changes in the structure of the 
vortex cores the vortices can actually be pinned at the minima of the potential. In Sec. II 
we perform numerical mean-field calculations, and find a sequence of first-order transitions 
between site-centered and plaquette-centered vortex lattices. In Sec. Ill we use a reduced 
basis ansatz to perform analytic calculations near the Mott boundary, and as a result show 
how the theory at the Mott boundary is related to the Hofstadter butterfly spectrum. In 
Sec. IV we summarize our results. 
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II. NUMERICAL CALCULATION OF VORTEX-LATTICE STATES 



A. Mean-field theory of the rotating Bose-Hubbard model 

We consider a deep lattice where we can make a tight-binding approximation jl4j |. and 
the system is described in the rotating frame by a Bose-Hubbard Hamiltonian jlol ]. 



Hrbh — —t 



{^d\cij exp 






i dr ■ A(r) 


+ H.c. j 









V^i (hi - 1) - fi^hi . (1) 
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Above, the operator d\ (a,) creates (destroys) a boson and rij is the number operator at 
optical-lattice site i. The subscript denotes a nearest-neighbor sum. The parameters t, 
U, and fi are the hopping matrix element, the on-site repulsion strength, and the chemical 
potential, respectively. Rotation produces the vector potential A(r) = (m/h) (ft x r) = 
ttu (xy — yx), where v is the number of circulation quanta (h/m, where m is the atomic 
mass, and h is 2ir times Planck's constant K) per optical-lattice site. Rotation also produces 
a harmonic centrifugal-potential which we have assumed is cancelled by a harmonic trap. 
Although we choose to work in the symmetric gauge our results are not gauge dependent. 
Scaling energies by U and distances by the lattice constant, the system is characterized by 
the unitless parameters t (= t/U), jl (= n/U) and v. 

We choose to model a uniform system, rather than explicitly considering a harmonic 
trap, because we feel that this approach gives more understanding of the phenomena. In 
addition, we also restrict ourselves to two dimensions, where the physics we are investigating 
is particularly clear. This geometry can be engineered by applying a sufficiently strong 
optical lattice in the z-direction which prevents hopping in that direction [if]]. Also, a 
rapidly-rotating BEC can assume a similar geometry through centrifugal distortion of its 
density profile [2]. Furthermore, we restrict ourselves to the case where the rotation speed 
is tuned so that v is a rational fraction, thus avoiding the commensurability issues which 
generically occur [k]]. In the strong optical-lattice limit, the vortex lattice will share the 



geometry of the optical lattice 
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As one approaches the superfluid-Mott boundary from weak coupling, the vortex cores 
evolve from empty to containing the Mott phase 15]. This happens because when the 
superfluid order is suppressed in the vortex core, the Mott phase is energetically favorable 
compared to the vacuum. This raises the possibility that the energy of the vortex lattice will 
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be reduced if the cores are centered on optical-lattice minima, "sites", rather than at the 
potential maxima, "plaquettes" . A competing effect is that if the vortices are site-centered 
then the overlap of atomic wavepackets centered at neighboring sites will be reduced, raising 
the kinetic energy. We find that the interplay between these effects leads to a rich structure. 

To model an infinite vortex-lattice we perform self-consistent Gutzwiller mean-field cal- 
culations on a two-dimensional square-lattice supercell made up of L sites per side, where 
each site is an optical-lattice potential minimum. We focus on the simplest case where 
each supercell contains one quantum of circulation, which produces a ground-state solu- 
tion containing one singly-quantized vortex per supercell, and v = (1/L 2 ). The Gutzwiller 
mean-field theory can be viewed as a variational calculation where one minimizes (Hrbh) 
over the Gutzwiller product-states [fjj], |^) = Y\ { (^2 n fn\ n )i)i where i is the site index, n 
is the particle number, and \n)i is the n-particle occupation-number state at site i. Mini- 
mizing (Hrbh) with respect to f % * with the constraint ^ n |/^| 2 — 1 = gives L 2 nonlinear 
eigenvalue equations, one for each site, 



-t ((ak}V^f m ^Rjk + (a{)Vm + lf m+1 R k ^ + (^m 2 - L + | j m + xA f j m = , 

k, nn of j \ \ / / 

(2) 

where the sum is over all nearest neighbors of site j, mis the particle- number index, Xj is a 
Lagrange multiplier, and Rjk = exp i dr ■ A(r) , where i = y/—l. We iteratively solve 
these equations: first choosing a trial order-parameter field |aj°^|, where atj = (a,j); then 

updating it by = J2 n V^fn-i ({ a i P ~^}) ^" (l^ - ^})' w ^ iere P ^ s ^ e iteration index. 
Similar calculations were performed by Oktel et. al. [17j (in a strip geometry) and Wu et. 
al. (in a square geometry) to produce vortex lattices. 

We perform calculations in the neighborhood of the n = 1 Mott phase, so the occupation- 
number distribution of each site will be peaked about 1, with small variance. Hence we only 
need to allow /J, f( and /| to be nonzero. In most cases f\ an d /o wm De much smaller than 
f(. We find that using a larger occupation- number basis causes slight shifts of the boundary 
curves and the energy differences between plaquette- and site-centered vortex-lattice states, 
but the position of the Mott-lobe is unchanged. To model the infinite vortex lattice with 
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Curve number 


2 


3 


4 


1 


0.034 


0.00087 


0.000017 


2 




0.012 


0.0013 


3 






0.0046 



TABLE I: Separation between each structural boundary curve (L = 1 — 4) and its corresponding 
Mott lobe, quantified by At at p, = v2 — 1 (the Mott-lobe tip, see Fig. [1]). Curve number 1 refers 
to the curve closest to the Mott lobe, curve number 2 is the next curve out, etc. 



our (LxL)-supercell we use magnetic boundary conditions 
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a(x + L,y) = a(x,y)exp -i-(2y -y) 

71 

a(x,y + L) = a(x,y)exp +i— (2x — x) 



(3) 
(4) 



where a (x,y) = (dj), and (x, y) are the Cartesian coordinates of site j, and (xo,2/o) are free 
parameters which correspond to the coordinates of the vortex core in our supercell. 

B. Results and discussion 



The phase diagrams for L = 1 — 4 are displayed in Fig. [TJ Each phase plot has the familiar 
obe-shaped Mott-insulator region in the deep- well limit, whose size varies as one changes n v 
I7I . Q. We refer to the plaquette-centered vortex-lattice phase by the symbol P, and the 
site-centered vortex-lattice phase by the symbol S. As shown in these phase diagrams we find 
alternating bands of P and S. Moving from weak (large t) to strong coupling (small t) we 
find for L = 1: P; L = 2: PS; L = 3: PSP; L = 4: PS PS; L = 5 (not pictured): PSPSP. 
The bands get very narrow as one increases L and as one approaches the Mott lobe. Table 
D] gives the width of the various phases along the line jl = jl c , where p, c (= a/2 — l) is the 
scaled chemical potential at the tip of the n = 1 Mott lobe. 

There are several important features of these phase diagrams. First, the outermost vortex- 
lattice phase is always P, since a shallow optical-lattice potential pins vortices to the maxima 
of the potential. Second, for the values of L we have investigated, the phase diagram of a 
square vortex-lattice configuration characterized by n„ = 1/L 2 has L phase boundaries. 
Third, the innermost vortex-lattice phase alternates between P (odd L values) and S (even 
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FIG. 1: (wide, color online) (a)-(b) Structural phase plots for the cases L = 1 and L = 2, re- 
spectively. Dimensionless parameters t = t/U and fi = ^/J/ represent hopping amplitude and 
chemical potential, respectively, where each quantity is normalized by the on-site interaction. The 
plot labels P, S and MI refer to P-centered, S-centered and Mott-insulating phases, respectively, 
(c) The L = 3 phase plot, where shading is used to emphasize the thin reentrant P phase, (d) A 
closeup of the critical region of the Mott lobe in (c); the reentrant phase is more clearly resolved, 
(e) The L = 4 phase plot, on this parameter range, the inner structural-boundary curve cannot be 
discerned from the Mott lobe, (f) A closeup of the critical region of the Mott lobe in (e); shading 
is used to resolve the second reentrant phase region (S phase). 

L values). And finally, the phase boundaries appear to share a universal hyperbolic shape. 
Although we have no explanation for the second observation, below we will explain the 
others. 

We analyze the nature of the vortex-configuration phase transition by studying how the 
energy depends on the location of the vortex core in a single supercell. Figure [2] illustrates 
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FIG. 2: (color online) Energy vs core placement. Vortex core position (xo,yo) i n units of optical- 
lattice spacing with (xo,yo) = (0,0) corresponding to a vortex centered on a site, and (xo,yo) = 
(0.5, 0.5) corresponding to a vortex centered on a plaquette. These plots correspond to the L = 3 
recurrent phase boundary at [i = \/2 — 1, and 0.0519 < t < 0.052. In (a) (t = 0.0519) and (b) 
(t = 0.052) the vertices of the red (gray) lines are sites, and the plots are shaded so that darker 
(lighter) corresponds to lower (higher) energy. Plot (a) [(b)] corresponds to the P (S) state for t 
just below (above) the boundary, (c) A composite of energy vs core-position curves on the diagonal 
line yo = xq £ (—0.5, 0) (from plaquette to site), for t between the spinodal points of the boundary 
For each curve E [xq] = [E (xq) — E (—0.5)] / Ey[ott, where E (xq) = {Hrbh) (xq). From top to 
bottom, this plot has 15 lines corresponding to t max = 0.051902 and t m j n = 0.0519015, with spacing 
At = 7.5 x 10~ 7 . 

that the transitions are discontinuous. We quantify the abruptness of the phase transition 
by measuring the width of the coexistence region; that is, we calculate the difference in t 
(at fixed p) between spinodal points where each of the two energy minima disappear. As 
shown in Table [Til the coexistence region becomes thinner as L increases, and as the system 
moves toward the insulating phase. 

The experimental consequences of our findings depend crucially on the energy difference 
of the two configurations. For example, the lattice will no longer be pinned if the temperature 
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L = 

Curve number 


2 


3 


4 


1 


0.014 


7.5 x 10~ 6 


2.8 x 10~ 7 


2 




0.0004 


1.5 x 10~ 6 


3 






2.5 x 10~ 5 



TABLE II: Coexistence region widths, A t, at ft, = y/2 — 1 (Mott-lobe tip) for the structural phase 
boundaries (L = 1 — 4). Widths are determined by finding the distance between spinodals. Curve 
number 1 refers to the boundary curve closest to the Mott lobe, curve number 2 is the next curve 
out, etc. 

T exceeds this energy. On the line jl = fi c we plot these energies in Fig. [31 The pinning 
energy decreases rapidly as L increases, and as the system approaches the insulating phase 
in parameter space. We find that the phases inside the outermost P phase all have tiny 
energy differences. To even see the L = 2 transition one requires a temperature below 
0.15 nK. Hence our findings are mainly of academic interest. If the temperature is large 
compared with the pinning energy then the vortex configuration will be determined by 
the competition between vortex-vortex interaction, which favors a triangular vortex-lattice 
phase, and entropy, which favors a disordered vortex-liquid. 

An additional concern is that the structures we find might be in part an artifact of the 
mean-field theory. Even if this is the case, we believe it is valuable to understand the struc- 
ture of the mean-field theory. Furthermore, in the following sections we will give arguments 
which suggest that those results are more general. Finally, we note that experiments are 
currently far from the regime we consider. 



III. ANALYTIC THEORY NEAR THE MOTT-BOUNDARY 

Very near the Mott phase we can linearize Eq. ([T]) and analytically calculate the state 
of the system. During preparation of this paper, Umucahlar and Oktel 2l|] presented an 
independent study with substantial overlap of this section. 



S 




FIG. 3: (color online) Energy difference between P and S states with respect to t at fixed 
values of L. (a)-(d) correspond to L=2-5, respectively. The dimensionless energy difference 
AE = (Ep — Es) /U, where Eptg) = {Hrbh)p(s)- The P-centered configuration is always favored 
in the outermost phase region. The energy differences decrease with decreasing n v (increasing L), 
and also as the system approaches the insulating region (decreasing t). These numbers suggest 
that, in practice, a homogeneous vortex-lattice configuration is unlikely to be thermally stable in 
any of the inner phase regions. 

A. Reduced-basis ansatz and Harper's equation 

It is simplest to illustrate this method by starting with the case of a uniform system 
which is not rotating (0 = 0). The expectation value of this Hamiltonian with respect to 
our Gutzwiller product state is 

(H)/N = -at\a\ 2 + l -{h 2 ) - (fi + (n) , (5) 

where a is the number of nearest neighbors, and N is the total number of sites. As one 
approaches the n-particle Mott lobe we can, as before, make the ansatz that the single-site 
wavefunction is of the form = f n -i\n - 1) + f n \n) + f n +x\n + 1), with (/ n _i, f n , f„.+i) = 
(ei, i/l — ef — e|, e^j , where e« are small. One can readily verify that the terms neglected 
in this ansatz are of higher order in e. Minimizing (H) one finds that the chemical potential 
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at which e,- becomes nonzero is 



t ± =(n-±(l + atfj±J \oH* - (n + at + \ . 



(6) 



In particular, the tip of the Mott Lobe is at p, c 
2n+ 1 + lyjn (n + 1) \ 
Adding rotation, the energy divided by U is 



y/n{n + \) — 1, at c 



(H RBH ) = -^2 {tij a * a j + c - c ) + Yl \ ~ + 2 ) ^0 



(7) 



where t ijf - = t exp mv f* 1 (xy — yx) ■ dr 



. Again, near the Mott lobe we write 



(/U fl fUi) = ^/l-\a t \H\\^ + \\\n X\ 



a, 



(8) 



where a = a+O (a 3 ), and X\ = (1 ~~ \f^X l ). Note that unlike our previous calculations, 
we do not need to restrict n v = v. 

Next we minimize with respect to A* to find 



(Hrbh) 



$5 a i a 3 + c - c ) + 1 1 ~ n 'i 



n — fx \ i2 T-i 



(9) 



where i^Mott is the energy-per-site of the n-particle Mott state, and we have neglected terms 
of order a 3 . Next we minimize with respect to at. In the case of the 2D square lattice we 



arrive at a symmetric-gauge Harper's equation 
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— a (x + l,y) exp [iiruy] — a (x — 1, y) exp [—mvy] — a (x, y + 1) exp [— mvx] 

— a (x, y — 1) exp [+i7rz/x] + ea (x, y) = , (10) 



where 

€ ~ in + 1 

A simple gauge transformation, otj = <Xj exp 



In — a I n — a 
1 — n 

1 + fJL 



iixv f n (xy + yx) ■ di 



(11) 



along with the as- 



sumption that a (x, y) = exp (i^)(3 (x), brings Eq. (fTOj) into the more familiar form 



p {x + 1) + (3 (x - 1) + 2 cos (2-kvx - 7) /3 (x) = e/3 (x) , 



(12) 
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FIG. 4: (wide, color online) The blue (light gray) surface in (a) is the mean-field Mott boundary 
of the Bose-Hubbard model at zero temperature for chemical potential fi = {0, 1}, and circulation- 
quanta per optical-lattice site v = {0, 1}. The red (dark gray) curve on this surface [and outlining 
the bottom edge of the spectrum in (b)] demonstrates how, at fixed fl (the value in the figure 
is jl = 0.2), the value of t is inversely related to the edge eigenvalues of the Hofstadter butterfly 
spectrum shown in (b). The black curve on the boundary surface [and in (c)] demonstrates how, 
at fixed v (in this case v = 1/4), the value of t is just a familiar Mott-lobe boundary in the 
(t, fx) -plane, as shown in (c). 

where the circulation density v = p/q is a rational fraction, and 7 is a wavevector set to 
ir/2q in Ref. 22j. The eigenvalue spectrum of Eq. ffTUl) has an intricate fractal-structure 
known as the Hofstadter butterfly 221 ] . 

Fixing v and pL, the corresponding point on the Mott lobe is the smallest t for which 
Eq. ( jTTT) is an eigenvalue of Eq. (|T0l) . This condition is satisifed by the largest eigenvalue 
e = e c dgc iy) of Eq. ( jTUl) . We call this largest eigenvalue the edge eigenvalue. The Mott 
boundary is then given by 

1 n — fi I n — fi 



t 



n- 



(13) 



CedgeM 1+1 \ 1 + \X, 

where n is the integer corresponding to the total-particle density in the Mott lobe. This 
remarkable relationship is illustrated in Fig. HI In the non-rotating case we find e e d ge [v = 0] = 
4, and Eq. ([TBI reduces to Eq. (jSJ). 

The eigenvectors of Harper's equation have rich topologies. The highest band (corre- 
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FIG. 5: (wide, color online) Hofstadter butterfly eigenvectors, for ^ = 1/100 in a 10x10 supercell. 
The position coordinates (xq, yo) are in units of the optical lattice spacing, and the order parameter 
density \a\ 2 is normalized so that over a single supercell ^(z j/o)l a ( x °' 2/o) 1 2 = 1- The bands are 
indexed with n = 1 for smallest central-eigenvalue, n = 2 for next smallest, etc. (a)-(c) Plots of 
order-parameter density |a| 2 for bands n = 100, n = 97 and n = 91 respectively, (d)-(f) The 
corresponding complex-phase fields. At each site is the base of an arrow pointing in the direction 
(Re [a] ,Im[a\), and with length proportional to \a\. Positively (negatively) charged vortices are 
labeled with a red "+" (blue "— "). The green boundary encloses one unit cell. The size and shape 
of this boundary are fixed, but varying e will shift its position. The n = 100 plot has a single vortex 
with charge +1. The n = 97 state has a central doubly-quantized vortex of charge 2, connected by 
domain walls to vortices of charge —1 near the faces of the cell. Vortices of charge +1 lie near the 
corners. The n = 91 pattern contains 8 "+" -vortices and 7 "—"-vortices in each unit cell. 
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sponding to the largest e's) contains states with regular arrays of singly quantized vortices. 
Changing e continuously changes the location of the vortices relative to the lattice. The 
lower bands include states with more complicated structures with multiple vortices of op- 
posite signs. Also, the band structure is symmetric with respect to reflection about e = 0. 
Illustrative structures are shown in Fig. 

The edge state corresponds to an array of singly quantized vortices. For v — 1/L 2 we find 
that for even or odd L these vortices are site-centered or plaquette-centered, respectively. 
This explains our previous observation of alternating vortex-lattice phases corresponding to 
even and odd L values. 



B. Discussion 



Why does the Hofstadter butterfly, a pattern associated with noninteracting particles, 
appear near the Mott lobe, where the interactions are strong? The answer is that near the 
Mott-lobe boundary most of the atoms are static, with only a dilute gas of mobile particles 
and/or holes. The diluteness of these excitations leads to single-particle physics. 

It should be noted that this explanation does not depend on the approximations of 
mean-field theory. Even including fluctuations, near the Mott lobe (with the exception of 
the region immediately about the tip), the system is described by a weakly- interacting gas of 



excitations 



20 



231 ] . Thus it is unlikely that the structural transitions we find are an artifact 



of mean-field theory. Interactions between the excitations can be included in our mean-field 



theory by including higher-order terms in equation 
one recovers a nonlinear Schrodinger equation 



If one approximates (nf) = |aij 



- £ 

j , nn of k 



i^otj exp 


". f 

ITXV / 







(xy — yx) ■ dr 



I ,o /i 



0. 



(14) 



We should mention that one can also study Hofstadter butterfly physics far from the 



Mott lobe by using a Feshbach resonance 



24] to tune the interaction of a gas of bosons 



trapped in a deep, rotating optical lattice. Merely reducing the lattice strength is probably 
insufficient, as the tight-binding approximation is apt to break down. 
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IV. SUMMARY 



We have analyzed vortex-lattice phases in a deep optical-lattice potential using the mean- 
field theory of the rotating Bose-Hubbard Model in a two-dimensional square-lattice at zero 
temperature. We observed several transitions between site-centered and plaquette-centered 
vortex states. For the (L x L)-supercell calculation (corresponding to n v = 1/L 2 ) there are L 
boundary curves - L — 1 structural curves, and the Mott lobe. We found that the structural 
transitions are discontinuous, and we quantify trends in the widths of the corresponding 
coexistence regions as well as trends in the spacing of the structural boundary lines in 
parameter space. The boundary curves share a universal hyperbolic shape. 

We also carried out an analytic study where we determined that the linear eigenvalue 
equation characterizing the Mott lobe also characterizes the Hofstadter butterfly spectrum. 
From this we determined an expression for the Mott-lobe boundary. This linearized analysis 
confirmed the vortex-core placement found in our numerical study. 
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